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I. INTRODUCTION 



The complex Langevin method solves in principle the 
sign problems arising in simulations of various systems, 
in particular in QCD with finite chemical potential. Af- 
ter being proposed in the early 1980s by Klauder [jjand 
Parisi 0, it enjoyed a certain limited popularity [3, H|, 
but very quickly certain problems were found. The first 
one was instability of the simulations with absence of 
convergence (runaways), the second one convergence to 
a wrong limit Nevertheless in recent years the 

method has been revived with sometimes impressive suc- 
cess • I n particular the use of adaptive stepsize has 
eliminated the problem of runaways [171 ]. 

But nagging problems remained due to the lack of clear 
criteria to decide when an apparently convergent simu- 
lation actually represented the truth. This was linked to 
the lack of a clear mathematical basis for the method, 
that would at the same time also provide criteria for its 
applicability. The purpose of the present paper is to clar- 
ify the situation at least to some extent. While we are still 
not able to close certain mathematical gaps and reach 
a complete analytic solution to the problems that have 
plagued the method, we give some strong numerical evi- 
dence that the method is correct in some cases and also 
suggest a plausible explanation for the failure in other 
cases; this leads to some pragmatic conclusions suggest- 
ing how to proceed in practice in a way that promises 
credible results. 

The paper is organized as follows. In Sec. [TTJ we give 
a formal justification of the method, highlighting the as- 
sumptions underlying the derivation. In Sec. IIIII three 
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main questions raised by the formal arguments are listed. 
We then focus on one particular issue, boundary effects, 
in Sec. IIV1 and present detailed case studies in Sec. |V] 
Tentative conclusions are given in Sec. IVI1 



II. FORMAL DERIVATIONS 



For simplicity we concentrate here on models in which 
the fields take values in flat manifolds A4 = K™ or 
M = T n , where T n is the n dimensional torus (S* 1 )" 
with coordinates (x\, . . . ,x n ). The complications that 
arise when the fields live in nontrivial manifolds, as is 
of course the case in QCD, have been successfully dealt 
with in the literature (see for instance Ref. [l8[ for real, 
Refs. pil [H, Gil H(| for complex Langevin dynamics). 
But these complications are not really relevant for our 
discussion. 

As is well known, the idea is to simulate a complex 
measure exp(— S)dx, with S a holomorphic function on 
a real manifold A4 , by setting up a stochastic process on 
the complexification Ai c of A4, such that the expecta- 
tion values of entire holomorphic observables O in this 
stochastic process converge to the ones with respect to 
the complex measure exp(— S)dx. 

The complex Langevin equation (CLE) on A4 C is 



dz = —VSdt + dw, 



(1) 



where dw denotes the increment of the Wiener process 
and the equation is to be interpreted as a real stochastic 
process, namely 



dx =K x dt + dw, 
dy =K y dt, 



(2) 
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with 



K x = -ReV x S(x + iy), 
K y = — liaS/ x S{x + iy). 



(3) 



A slight generalization of Eq. ^ that has been consid- 
ered and will play a role in this investigation is 



dx =K x dt + y N R dw R , 
dy =K y dt + \/ Njdwi, 



(4) 



where dw/j and dui/ are independent Wiener processes, 
iVj > and = Nj + 1. This is usually referred 
to as complex noise. The introduction of a nonzero iVj 
makes it possible to solve the Fokker-Planck equation 
(see below) numerically and also allows a random walk 
discretization of the complex Langevin process: 



Sx(t) = ±U) X , P x> ± = i ( I ± t an h 
<fy(t) = ± uiy, P y> ± = - I 1 ± tanh 



2JV fl 



3/' ± a, j- 2 

W x = y/2N R 6t, L0 y = y/2N!St, 



^Ky 



(5) 

(6) 
(7) 



where P± are the transition probabilities and we have 
defined the steps ui such as to have the same St in both 
sub-processes, to ensure correct evolution. 

By ltd calculus, if / is a twice differentiable function 
on M c and 



z(t) = x(t) + iy(t) 



(8) 



is a solution of the complex Langevin equation ((4|), we 
have 

±(f(x(t),y(t))) = (Lf(x(t),y(t))), (9) 

where L is the Langevin operator 

L = [N R W X + K x ] V x + [NiVy + K y ] V v , (10) 

and (/) denotes the noise average of / corresponding to 
the stochastic process described by Eq. ((H). In the stan- 
dard way Eq. Q leads to its dual Fokker-Planck equa- 
tion (FPE) for the evolution of the probability density 
P(x,y;t), 



d_ 

di 



P(x,y;t) = L T P(x,y;t), 



with 



L T = V x [N R V X - K x ] + V v [AW, 



(11) 



(12) 



L T is the formal adjoint (transpose) of L with respect to 
the bilinear (not hermitian!) pairing 



(f,P) = / f(x,y)P(x,y)dxdy, 



(13) 



(Lf,P) = (f,L T P). (14) 
Note that the FPE has the form of a continuity equation 



8 . 

gj.P{ X 1 J/5 *) = ^xJx + Vy Jy, 



(15) 



where (J x , J y ) is the probability current in the 2n dimen- 
sional space M c , given by 

J x = (JVkV, - K x ) P, Jy = {NiV v -K v )P. (16) 

We will also consider the evolution of a complex density 
p(x) on A4 under the following complex FPE 



d 

— p{x:t) = L%p(x;t), 



(17) 



where now the complex Fokker-Planck operator Lq is 



Ll = V, [V, + {V x S(x))} 



(18) 



A slight generalization will be useful: For any yo 6 A4 
we consider a complex Fokker-Planck operator Ly given 
by 



L yo= V * [V x + (V x S(x + iy ))} 



Ly o is the formal adjoint of 



L vo = [V? - (VxS(x + iy ))} V x . 



(19) 



(20) 



The operators Ly act on suitable complex valued distri- 
butions (measures) on M. , parameterized by the real vari- 
ables (xi, • • • , x n ). But they do not allow a probabilistic 
interpretation, because they do not preserve positivity. 

For any y Eq. (fT7| with Lq replaced by Ly g has the 
complex density 

p Vo (x;oo) <xexp[-S(x + iy )] (21) 

as its (hopefully unique) stationary solution. 

We next consider expectation values. Let O be an 
entire holomorphic observable with at most exponential 
growth; then we set 

//m _ JO(x + iy)P(x,y;t)dxdy 



J P(x,y;t)dxdy 



and 



_ / Q(x)p(x; t)dx 
{U)pW = Jp(x;t)dx ■ 

We would like to show that 

(0) m = (0) p(t) , 

provided the initial conditions agree, 

(C)P(O) = (C)p(O); 



(23) 

(24) 
(25) 
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which is true if we choose 



P(x,y;0)=p(x;0)5(y-y ) 



(26) 



(for any yo). In the limit t — > oo the dependence on the 
initial condition should of course disappear by ergodicity. 

The goal is to establish a connection between the 'ex- 
pectation values' with respect to p and P for the class 
of observables chosen (entire holomorphic with at most 
exponential growth). The idea is to move the time evo- 
lution from the densities to the observables and make 
use of the Cauchy-Riemann (CR) equations. Formally 
(i.e. without worrying about boundary terms and exis- 
tence questions) this works as follows: first we use the 
fact that we want to apply the complex FP operators 
L yo only to functions that have analytic continuations to 
all of M c - On those analytic continuations we may act 
with the Langevin operator 



L= [V 2 -(V*S(z))]V 2 



(27) 



whose action on holomorphic functions agrees with that 
of L, since on such functions V y = iV x and A x = — A y 
so that the difference L — L vanishes. 

The proliferation of Langevin/Fokker-Planck operators 
may be somewhat bewildering, but it is important to 
realize that L,L,L yo are really all different operators: 
while L and L act on functions on A4 C (i.e. functions of 
(xi,yi, . . . ,x n ,y n )), agreeing on holomorphic functions, 
but disagreeing on general functions, L yo acts on func- 
tions on M., i.e. functions of {x\, . . . , x n ). 

We now use L to evolve the observables by the equation 

d t O{z; t) = LO{z- 1) (t > 0) (28) 

with the initial condition O(z;0) = O(z), which is for- 
mally solved by 



0(z;t) = exp[tL](D(z). 



(29) 



In Eqs. ([281 [29]) . because of the CR equations, the tilde 
may be dropped, and we will do so now. So we will have 



d t O(z;t) =LO{z;t) (t > 0), 
with its formal solution 



(30) 



In fact Eq. 
tions 



0{z;t) = exp[tL]<D(z). (31) 
is also equivalent to the family of equa- 



d t O(x + iyo;t)=Ly 0(x + iy ;t) (t > 0) Vy - (32) 

The hrst thing to notice is that 0(z; t) will be holomor- 
phic if 0(z; 0) is. This can be seen as follows: Let dj 
(i = j, . . . , n) be the Cauchy-Riemann operators defined 
by 



dj = d x 



id. 



Mi ■ 



(33) 



Applying this to both sides of Eq. (l28l) and using the 
fact that by the holomorphy of S dj commutes with L, 
we find 



d t djO(z; t) = Ld 3 0(z; t) (t > 0) ; 



(34) 



this is just Eq. (|2"51) again with 0(z;t) replaced by 
BjO(z; t). Under the assumption that L generates a semi- 
group acting on BjO(z; t), Eq. _(l34l) has a unique solution; 
since the initial condition is djO(z;0) = we conclude 
that BjO(z,t) = for all t > 0, i.e. 0(z;t) satisfies the 
CR equations for alH > and all j = 1, ... , n. So 0(z; t) 
is holomorphic in each component Zj separately. By Har- 
togs' theorem [U [22j this implies joint holomorphy. 



We now consider, for < r < t, 



F(t,r) = J P(x, y;t- t)0(x + iy; r)dxdy, (35) 

and claim that it interpolates between the p and the P 
expectations: 

F(t,0) = (0) P(t) , F{t,t) = {0) p{t) . (36) 

The first equality is obvious, while the second one can be 
seen as follows, using Eqs. (f2"6l l3"Tj) . 



F(t,t) = / P(x,y;0) exp(tL)(D{x + iy;0)dxdy 



= J p{x;0) (exp(tL yo )G(x + iy ;Q))dx 
= J(D{x + iy Q ; 0) (exp(tL£>(x; 0)) dx 



(37) 



where we only had to assume that we can integrate by 
parts in x without worrying about boundary terms. 

Our desired result follows if we can show that F(t, r) 
is independent of t. To sec this, we differentiate 

^-F{t, t)=- f (L T P{x, y; t - r)) 0(x + iy; r)dxdy 

+ / P(x,y;t - t)LO(x + iy;r)dxdy. (38) 



Integration by parts then shows that the two terms can- 
cel, hence -^F(t, t) = and thus 



p) P{t) = (0} p(t 



(39) 



It is important to notice that this holds for all Nj; 
whereas the left hand side seems to depend on Nj, the 
right hand side is manifestly independent of it. 

If we knew in addition that 



lim(0) p(t) = (0)„ (oo) , 



(40) 



with p(oo) given by Eq. (|21[) with yo — 0, we could now 
conclude that the expectation values of the Langevin pro- 
cess relax to the desired values; this convergence would 
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follow if we knew that the spectrum of Ly o lies in a half 
plane Re z < and is a nondegenerate eigenvalue. But 
note that we do not really need convergence of P(x, y; t) 
for Eq. (|4"P|) to hold, since it will only be tested against 
analytic observables O. 

Nevertheless the numerical evidence in many cases 
points to the existence of a unique stationary probability 
density P(x,y;oo). The corresponding probability cur- 
rents are divergenceless, but unlike the situation in the 
real Langevin process, they cannot vanish. A general fea- 
ture of the stationary distribution that can be read off 
the FPE is the following: Assume that (x , yo) is a local 
stationary point of P, then 

(N R A X + N!A y )P = V X {K X P) + Vy(K y P), (41) 

for x — XQ,y = yo. So if Ni > 0, a local maximum 
of P can only occur where the divergence of the drift 
force is negative and a minimum where it is positive. 
For Nj — the conclusion is even stronger: where the 
divergence is negative (positive), there cannot be a local 
minimum (maximum) in x for fixed y. These properties 
provide some checks on numerical solutions. 



III. QUESTIONS 

There are three main questions raised by the formal 
arguments in the previous section: 

(1) Can the operators L, L, L ya and their transposes be 
exponentiated; in more mathematical language: do these 
operators generate semigroups on some suitable space of 
functions? 

(2) Are the various integrations by parts justified, 
which underlie the shifting of the time evolution from 
the measure to the observables and back, or are there 
boundary terms to worry about? 

(3) Are the spectra of L,L yo and their transposes en- 
tirely in the left half plane and is a nondegenerate eigen- 
value? 

Concerning the fir st q uestion, there are treatises (see 
for instance Refs. [23|,[24]]) giving rather general sufficient 
conditions for the existence of a semigroup generated by 
differential operators of the general type considered here. 
Unfortunately it seems that the cases we have to deal 
with here are not covered by those general results; the 
main difficulties are (1) the strong growth of the drift 
given by the gradients of the action in some complex 
directions and (2) the fact that the drift is not always 
'restoring'. 

Question (1) for L, L T is intimately related to the ques- 
tion whether the stochastic process given by the complex 
Langevin equation exists for arbitrary long times. This is 
not obvious because typically in the classical (no noise) 
limit there are trajectories that go to infinity in finite 



time. While those trajectories occur only for a subset of 
measure zero of initial conditions, it is not obvious what 
happens after adding the noise. On the one hand, the 
noise will typically kick the process away from the un- 
stable trajectories; on the other hand it may also kick 
it near the unstable trajectory, inducing very large ex- 
cursions of the process. We will later illustrate this with 
some examples. 

But let us say that the accumulated numerical evidence 
points not only to the existence of the process for arbi- 
trarily large times, but also to the existence of a unique 
equilibrium measure for the process; unfortunately we 
could neither find results in the mathematical literature 
that would imply this, nor could we prove ourselves that 
this is the case. 

Concerning the exponentiation of L, L yo , which should 
be easier, we still could not establish mathematically that 
it is possible on a space of functions containing the most 
obvious observables, such as exponentials. 

There is a useful criterion for the existence of a 
bounded semigroup generated by an operator A on a 
Hilbert space [231 ] : A generates a bounded semigroup if 
it is dissipative, i.e. if A + A> < 0. Unfortunately even 
in the simplest cases of a quadratic S the corresponding 
Langevin and FP operators are not dissipative. So they 
can at best generate exponentially bounded semigroups; 
if in addition the spectrum is in the left half plane, con- 
vergence to the equilibrium should still take place. 

For the second question it would be necessary to have 
good control over the falloff of the solutions of the FPE 
in the imaginary directions: if we insert the observ- 
ables Ok(z) = exp(ifcz) into Eq. (|39|) . we get for the 
Fourier transform (Fourier coefficients in the compact 
case) p(fc; t) of the complex density p(x; t) 



m t) = j^y t J P(x, y; ty k *- k ydxdy. (42) 



This makes sense for all k only if P(x, y; t) decays more 
strongly than any exponential in imaginary direction. 
Our case studies described in the following sections in- 
dicate that this does not seem to be the case: in our 
first example the decay is probably exponential, but not 
stronger; in our second example the decay seems to be 
even weaker (0(\y\~ r ) with r sa 2), so that exponentials 
cannot be used as observables. The remainder of the pa- 
per is mainly devoted to studying Question (2) in some 
toy models. 

Finally let us remark that the third question is more 
difficult than the first one, and again the answer is not 
known rigorously. But again the numerical evidence 
strongly suggests a positive answer, depending on the 
model and the parameter values, in many interesting and 
relevant cases (see e.g. Refs. (HQ). 
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IV. BOUNDARY EFFECTS 



In this paper we are mainly concerned with Question 
(2). Even though we have very little analytic control, 
careful numerical studies reveal that, as remarked, the 
answer is generally no! In our case studies we find in- 
dications that the probability density P(x,y;t) indeed 
relaxes to an equilibrium density P(x,y;oo), but that 
that limiting density decays at best like an exponential 
for \y\ —¥ 00, in other cases only power-like. This lim- 
its the class of observables for which the integrations by 
parts can be performed without boundary terms. 

But let us first take a closer look at the integrations 
by parts that occur in the formal arguments of Sec. |n] 
The danger lies in a possibly insufficient falloff in the 
imaginary (y) directions, whereas in the real (x) direction 
we have either compactness or sufficient falloff due to the 
behaviour of the action S. 

Let us remark that for Nj > the operators L and L T 
are uniformly strictly elliptic; this is important, because 
it implies regularity for any solution of the stationary 
FPE 25]. It is also to be expected that the semigroup 
exp(tL T ) has a smooth (even real analytic) kernel so that 
for t > P(x,y;t) will be smooth (real analytic). This 
is supported by our numerical studies. So the problem- 
atic point is to show that the two terms on right-hand 
side of Eq. (|38|) actually cancel. This may fail because 
the observables O typically grow in imaginary direction, 
whereas the decay of P(x, y; t) (always assuming it exists) 
may be insufficient to compensate for it. 

Let us see in a little more detail how the argument 
for the independence of Nj may fail. For simplicity of 
presentation we consider the one-dimensional case (n = 
1). We write the Langevin operator pH|) as 



L = L Nl=0 + NjA, 



(43) 



where Ln i= q is the Langevin operator for Nj = 0, Nr 
1. Then let us consider 



d 



8 



dm 

We use the formula 



P(x, y\ t)0(x + iy)dxdy 



{e t(L "i= 0+NlA) P{x, y; 0))O(x + iy)dxdy. 



(44) 



d 



= f dre TLT Ae^ LT , (45) 
Jo 

and integration by parts to rewrite Eq. (|4"4")l as 

J dr J P{x,y;t - t)AO(x + iy;r)dxdy + X. (46) 

The term denoted by X collects possible boundary terms 
arising in the integration by parts. It vanishes only if the 



decay of P(x, y;t—r) is strong enough to offset any possi- 
ble growth of 0{x, y; r); otherwise it may either converge 
to a finite nonzero value or diverge. By the CR equations 
the first term vanishes, but the uncontrolled boundary 
term X remains. 

Let us look at a simple example that shows how and 
when the formal argument fails. Consider 

* = /<An«;<))0(*+«^ 



P(x, y; t)AO{x + iy)dxdy 



(47) 



(where the second term vanishes on account of the CR 
equations). By the formal argument this would be zero, 
being just a boundary term, but careful application 
of integration by parts, at first over the finite domain 
— F_ < y <Y + , gives 



X = lim dx (y y P(x,y;t))0{x + iy) 

Y± — ► oo 



-P{x,y;t)V y O{x + iy) 



-y_ 



(48) 



Here it is clear that what matters is the combined asymp- 
totic behaviour of P and O, and depending on the ob- 
servable, X may be zero, finite, or divergent. Of course 
the form of the boundary terms is less simple when using 
integration by parts in Eqs. (|38I44I) . but we expect that it 
is still the decay of the products like PO, OV y P, PV y O 
that is relevant. 

When trying to investigate the effects of the boundary 
numerically, one would in principle like to use the proba- 
bility density obtained without a cutoff. In practice this 
is, however, not feasible, and we therefore introduce a 
cutoff in the imaginary direction, which is not sent to 
infinity. Such a device is necessary for the solution of 
the FPE, and even though the CLE does not require it, 
for the purpose of comparison we also introduce it there. 
This will, however, introduce additional problems with 
the formal arguments relying on the CR equations as 
well as integration by parts. 

Concretely we proceed as follows: We restrict each yj 
to lie between — Y_ and Y+ and impose periodic bound- 
ary conditions on both the observables and the probabil- 
ity densities. This has a number of consequences. Firstly, 
observables 0(x + iy) will in general not be continu- 
ous across the 'seam', where we identify yj = —Y- with 
yj = Y + . They can therefore not be interpreted as con- 
tinuous functions and a priori the Ito formula © does not 
hold (it may still hold in the sense of distributions, which 
should be sufficient for our purposes). Furthermore, the 
jump across the seam will mean that the CR equations 
are no longer satisfied everywhere. For the evolved ob- 
servables the CR equations cannot be expected to hold 
anywhere exactly, as the violation that occurred initially 
only at the boundary gets propagated everywhere by the 
Langevin evolution. Similarly, the drift in the FPE is 
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expected to be discontinuous at the seam. Therefore we 
have to expect that P has a jump there as well; this in 
turn forces us to interpret the FPE in the sense of distri- 
butions. 

One might wonder whether it would not be better to 
limit the fluctuations in the imaginary direction by intro- 
ducing a smooth cutoff; but since such a smooth cutoff 
function will necessarily be nonholomorphic it will de- 
stroy the formal arguments even more; we therefore stick 
with the simplest choice of a periodic cutoff. 

We conclude therefore that the introduction of a cutoff 
— Y- < yj < Y + and imposing periodic boundary condi- 
tions leads to a breakdown of the formal arguments given 
in Sec. [Ill Although it is difficult to quantify precisely 
the effect of this, it seems reasonable to expect that it is 
still the behavior of PO and similar products at large \yj\ 
that determines what is happening. For the CLE it is also 
clear that a very large cutoff will practically not be felt, 
because the system very rarely will make contact with it. 
This is borne out by our numerics which clearly shows 
convergence to the limit of infinite cutoff. For the FPE, 
on the other hand the issue is less clear, because there 
are very large boundary terms arising from the gradients 
of the drift across the 'seam'. In any case, for the FPE 
we cannot directly compare with the cutoff-free results, 
because these do not exist. 



V. CASE STUDIES 

A. The U(l) one-link model 

To understand in more detail how boundary terms af- 
fect the behaviour of complex Langevin simulations, we 
studied in some detail the U(l) one-link model in the hop- 
ping (HDM) approximation that was already discussed in 
Ref. [Ill for Nj = 0. 

The action is 

S = — f3 cos z — kcos(z — iju) = — acos(,2 — ic), (49) 
with 



a = \J{P + kp>)(/3 + Ker^O, 

1, P + ne" 

c = t: hi t; • 

2 (3 + ne-v 

The complex drift force is correspondingly 

K = —S' = — (3 sin z — k sm(z — iji) 
— — a sin(z — ic), 

and the two components of the drift read 

K x = —Re S' = —a sin x cosh(y — c) , 
K y = — ImS 1 ' = —acosxs\nh{y — c). 



(50) 



(51) 



(52) 

(53) 
(54) 



As discussed in Ref. [12J there are two fixed points at 
(x, y) = (0, c) and (x,y) = (tt,c); the first one is attrac- 
tive, the second one repulsive. 

A special feature of this model is that for y = c the 
drift is purely in x direction. If in addition Ni = 0, 
the Langevin process will never leave the line y = c if it 
starts there (we emphasize that the properties discussed 
in this paragraph do not hold for the full U(l) one-link 
model, which was studied in detail in Ref. [12j). It is 
therefore straightforward to find an explicit solution to 
the stationary FPE, 



P(x,y; oo) oc e' s{x+lc) 5{y - c). 



(55) 



It follows that this model is actually equivalent to one 
with a real action, once we shift y — > y + c and replace /3 
by a (a and c now embody the dependence on the param- 
eters of the model). The numerics presented below show 
that the line y = c is an attractor for the Langevin pro- 
cess; this indicates that the solution (f55|) is unique (with 
proper normalization). These properties imply that for 
Nj = the dynamics is completely understood. 

When Nj > 0, the presence of the repulsive fixed 
point is responsible for the occurrence of large excur- 
sions, which are well known to be the scourge of complex 
Langevin simulations. For large y the drift terms dom- 
inate over the noise, and the Langevin process is essen- 
tially just a deterministic motion; the 'classical' trajec- 
tories are given by 



z(t) = i In 



1 - iCe~ 
1 + iCe~ 



(56) 



where the complex integration constant C is related to 
the starting point by 



C = tan 



2(0) 



(57) 



It is easy to see that all trajectories, except those starting 
on the unstable trajectories (lmz(0) = 7r ) are attracted 
to the stable fixed point at z = ic. 

As an illustration we show in Fig. [T] a scatter plot of 
a Langevin simulation clearly exhibiting the classical or- 
bits. There were 500 update steps between consecutive 
points and we used Nj — 1. To enhance the classical 
features, this simulation was done with the noise terms 
in the CLE (see Eqs. (|59l |60|) below) suppressed by a 
factor 10; this is of course equivalent to replacing e by 
e/100 while multiplying the force terms by a factor of 
100. This scatter plot should capture some generic fea- 
tures that will also be present for other choices of the 
model parameters. 

In order to numerically study the role of boundary 
terms and large \y\ values at nonzero Ni, we introduce 
a cutoff in imaginary direction, placed symmetrically 
around y = c, i.e. —Y + c < y < Y + c, and impose 
periodic boundary conditions. To see the effect of this 
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FIG. 1: Scatter plot for the U(l) one-link model at /3 = 1, 
K = 0, iVj = 1 with reduced noise (see text). 



cutoff, we compute numerically the expectation values 
of the observables exp(ifcz) for k = ±1,±2, for various 
values of Nj and Y, both by simulating the Langevin 
process and by numerical solution of the FPE. Note that 
these observables grow exponentially at large y. The ex- 
act values are given by 



I ikx\ _ I k{a) kc 
{ } ~Io(a) e 



(58) 



where Ik (a) are the modified Bessel functions of the first 
kind. 

The Langevin process is discretized in the usual way, 



) + V^Vfi^.n (59) 
y n+1 = y n + eK y (x n , y n ) + \J eNi rj Vtn (60) 

where r] XlV , and r\ y ^ ri are pseudorandom numbers with zero 
mean and variance 2. We use periodic boundary condi- 
tions in the imaginary direction, as stated above. We 
also use an adaptive step size [13], choosing e such that 
the product of e\K x + iK y \ < 0.001. To estimate the sta- 
tistical error we run 100 trajectories with independent 
random starting points. 

To solve the FPE numerically, we employ the period- 
icity in x and consider the Fourier decompositions 



Zn 



P(k,y;t) 



dx 
2^ 



e^P(x,y;t), 



(61) 
(62) 
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FIG. 2: Distributions P(x,y;t — > oo) in the U(l) one-link 
model, obtained from a numerical solution of the real FPE, 
for various values of Nr. Ni = 0.0001 (top), 0.01 (middle), 
0.1 (bottom). See the main text for further details. 



with the inverse transformations given by 

oo 

p(x;t)= J2 e- ik *p(k;t), (63) 

k— — oo 
oo 

P(x,y;t) = J2 e- tkx P(k,y;t). (64) 



8 




G-oN=0.0001 
H-H Nj=0.01 
0-ON=O.O5 
/WiN=0.1 



FIG. 3: Ni dependence of Re(e lz ) (lower points) and 
Re (e~ lz ) (higher points) from FPE for various values of the 
cutoff Y. The bottom figure zooms in on smaller values of 
Ni. The lines are guides to the eye, the horizontal dotted 
lines indicate the correct results (0.483564 and 0.592966 re- 
spectively). 



The FPE can be rewritten in terms of these modes as 



dtP(k, y; t) = (~N R k 2 + Njd 2 y )P(k, y- 1) 
+ | cosh(y - c) [(fc + l)P(k - 1, y; t) 



(k-l)P{k + l,y;t) 



+ - sinh(t/ - c)d y [P(k -l,y;t) + P(k + l,y; t)\ . (65) 

This equation is solved numerically with a discretized 
time step e = 10 -5 and a spatial discretization in y of 
S = y/e. We vary the cutoff Y in the y direction from 10<5 
up to BOOS (15005 in some cases). In the x direction, with 
— 7r < x < ir, we use 20 — 30 points. We found that con- 
vergence was reached after 10 6 time steps, corresponding 
to a Langevin time t ~ 10. We found convergence for all 
the values of Nj and Y studied. For Nj — 0, we were 
not able to solve the FPE numerically, due to the singu- 
lar behaviour, but the solution is known analytically, see 




0.5 0.75 

tanh(Y/1.581) 



FIG. 4: Cutoff (Y) dependence of Re(e iz ) (lower data) and 
Re (e _lz ) (upper data) for various values of Ni, for FPE (open 
symbols) and CLE (full symbols, note that the errorbars are 
much smaller that the points). The lines are guides to the 
eye. The horizontal dotted lines indicate the correct results. 



Eq. p)J). 

We fix the parameters of the model to be 
,3=1.0, k = 0.25, fi = 0.5. 



(66) 



In Fig. [5] we show some examples of the real probability 
distribution P(x,y,t) at large Langevin time t ~ 20 for 
various Ni and a given cutoff Y — 0.474. The y direction 
goes from left to right and the compact x direction from 
back to front. We observe that at small Nj = 0.0001, the 
distribution resembles the analytic result at Nj = 0. The 
distribution is very narrow in the y direction and bound- 
ary effects are not expected to play a role. Increasing 
Nj results in a wider distribution, and boundary effects 
become clearly visible. The apparent non-smooth be- 
haviour at the edge \y — c| = Y is to be expected from 
the discontinuity across the seam, discussed at the end 
of the previous section. 

After obtaining the distribution, expectation values of 
observables follow from Eq. (|2"2"|) . In Tables 1-7 (see 
end of paper) we compare the results of the Langevin 
simulation and the FPE for increasing values of Ni from 
up to Nj = 0.1. Note that the Langevin equation can 
be solved without a cutoff (Y — oo). The imaginary 
parts of the observables are consistent with zero. The 
data of the tables are also summarized in Figs. El H] (the 
FPE and the CLE data are indistinguishable at the scale 
of the figures; the rightmost points in Fig. 2] correspond 
to Y = 4.74 for FPE and oo for CLE). 

The following facts can be inferred from these results: 

(1) CLE and FPE give rather similar results, but 
sometimes they differ by several a (statistical error of 
the CLE simulation). 

(2) All the data show a clear dependence on Nj, in 
contrast to the conclusion of the formal arguments. 



9 



For larger Nj values both CLE and FPE give results 
different from the exact values. 

(3) The best results are generally obtained for the 
smallest 2Vj. In this case there is also the weakest Y 
dependence, in fact no Y dependence whatsoever for 
Ni = 0. 

Obviously the presence of the cutoff and periodic 
boundary conditions affects the CLE and FPE in a simi- 
lar way. But the Nj dependence shows the failure of the 
formal argument, even for Y — oo. At least for k = ±1 
one has a clear case of 'convergence to a wrong limit'. 
The data of the CLE with Y = oo and k = ±2 are actu- 
ally not really converged, except for Nj = 0. Observing 
the scatter plots for different Langevin times it appears 
unclear whether the observables exp(±2iz) really reach 
an equilibrium distribution - they tend to drift out to 
infinity, whereas their averages, while remaining small, 
suffer from huge fluctuations. On the other hand the ob- 
servables exp(±iz) seem to reach a stable distribution for 
Y = oo and any A?j. 




X 



FIG. 5: Scatter plot in the GP model with = 1, Ni = 1. 



B. The model of Guralnik and Pehlevan 



Guralnik and Pehlevan [15f studied an instructive toy 
model on Ai = R and Ai c — C, called GP model hence- 
forth. Its action is 



S' = - ;.;(z + iz 3 



(67) 



and was studied in connection with PT invariant but non- 
hermitian Hamiltionians (where PT indicates the com- 
bined action of parity and time reversal) [26j. 

The action (l67l) leads to the drift forces 



Kr 



-20xy, K y = 0(1 + x 2 - y 2 ) 



(68) 



There is a stable fixed point at z = i and an unstable 
one at z = —i. The 'classical trajectories' obtained by 
leaving out the noise are given by 



*(*) = 



z + i tanh(/3£) 
1 — izo tanh(/3t) 



(69) 



Since this is a Mobius transformation from w = tanh(/3i) 
to z(t) the trajectories are circles. They can be imag- 
ined to emerge from the unstable fixed point z = — i at 
t = — oo and go to the stable fixed point z — i as t — > oo. 
Those classical trajectories again can be seen clearly in 
the large excursions, since there the noise becomes negli- 
gible. Fig.[S]shows the result from a Langevin simulation 
at Ni = 1.0 and = 1.0 (in this case there are 50000 up- 
date steps between two consecutive points). 

Guralnik and Pehlevan give the exact results of the 



1.18 



0.2 0.4 0.6 0.1 

N, 



FIG. 6: Ni dependence of Im (z) in the GP model at = 1. 
first three moments at = 1, they are 

/-v Am) 



Ai(l) 



= -i, 

(z 3 )=i-( 



1.17631 



-0.1763z. 



(70) 



They solved the discretized Langevin equation numeri- 
cally, using Ni — 0.001, and obtained good agreement 
with those exact results. We did some more and prob- 
ably longer simulations at Ni = 0, 0.1, 0.5, 1.0. The 
results for Im (z) are shown in Fig. [6J again they show a 
clear dependence on iVj, in conflict with the formal rea- 
soning. While for small Ni there is agreement with the 
exact result, for larger Ni we again have convergence to 
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FIG. 7: Scatter plot in the GP model for /3 = 1, Nj = 0. 



the wrong limit. 

The results for Re {z 2 } show a similar behaviour, ex- 
cept that at iV> = 1.0 the result is -0.951(44), with a 
statistical error that is huge compared to the error 0.0007 
found at Ni — 0.5. The data for Im (z 3 ) show an even 
more dramatic failure for larger Nj: they diverge for 
Nj = 1.0, whereas for Ni = 0.5 fluctuations are large. 
Finally we measured Im(e lz ), which has an exact value 
of 2.2624. Here divergence becomes manifest already for 
Ni = 0.5 (but it might occur for all Ni > 0). 

In order to have an idea of the equilibrium distribu- 
tion for Ni = we show in Fig. [7] a scatter plot of 50000 
configurations in the complex plane. Non-Gaussian be- 
haviour is quite clear from this plot. Noticeable is the 
appearance of sharp edges of the distribution, possibly 
indicating jumps, but no 5 functions like in the U(l) 
case in the hopping expansion. To obtain these pictures 
we sampled over 60000 points taken at equal intervals of 
0.5 in Langevin time. Similar distributions have been ob- 
served in Ref. [l2j for the full U(l) one-link model. The 
non-Gaussian character is further demonstrated in Fig. 
[SJ the histograms for both Re z and Im z deviate strongly 
from a Gaussian distribution. 



VI. TENTATIVE CONCLUSIONS 



A. Explanation of wrong results 



We will now try to interpret these findings. They are 
quite analogous in the two model cases and we will try to 
reach some conclusions that can be generalized to more 
realistic models. 



1000 



:.' 000 



000 




FIG. 8: Distribution of Re z (top) and Im z (bottom) in the 
GP model with /3 = 1, Ni = 0. 



The question of convergence vs. divergence apparently 
depends on the values of iVj, but it is more plausible 
that there is no qualitative difference between the differ- 
ent positive values of Ni , only the time needed to observe 
the asymptotic behaviour is different. On the other hand 
for Ni = there seems to be really a qualitative differ- 
ence: the distributions develop discontinuities or even 8 
functions; but more important is the fact that they seem 
to drop very rapidly in the imaginary direction. 

We tentatively conclude that for iVj = the systems 
relax to equilibrium measures that show at least expo- 
nential decay in imaginary direction (in our simple U(X) 
model this decay is of course much stronger - the measure 
is zero for \y\ > c). 

For Ni > the situation is less clear, but it seems that 
in the U(l) model we get a decay at least like exp(— \y\), 
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but probably not stronger than any exponential. In 
the model of Guralnik and Pehlevan the data suggest 
a power- like decay (the power appears to be near 2). 

To sum up our tentative conclusions: 

There is a unique equilibrium distribution P(x,y; oo) 
for the Langevin processes for any Nj > 0, but for Nj > 
it shows limited decay, especially in the imaginary direc- 
tion. 

The type of decay for Nj > depends on the model 
considered; but for some observables their growth in 
imaginary direction may conspire with the falloff of the 
equilibrium measure in such a way that we obtain con- 
vergence to a wrong limit. In any case one should not 
expect convergence of the mean value for all holomorphic 
observables. 

For very small Nj and limited simulation time, the 
behaviour is of course indistinguishable from the one at 
Nj = 0, and one may reach a quasi-convergence to the 
right limit, even if an infinitely long simulation would 
diverge. 

If we try to generalize boldly from our toy U(l) model 
to real lattice gauge theories, we expect that for Nj > 
and for most interesting observables, such as Wilson 
loops, Polyakov loops etc., we have to expect boundary 
terms contributing even as the boundary is sent to infin- 
ity. For multiply charged loops and Nj > the situation 
could be even worse: those boundary contributions may 
diverge as the boundary is moved to infinity. But these 
problems probably will not occur for Nj — and do not 
show up at very small values of Ni either, at least if 
simulations are not run excessively long. 

We realize that definite conclusions about the falloff of 



the probability density P(x,y) in the y direction have not 
yet been reached; we intend to return to a more detailed 
study of this question in a future paper. 



B. Practical conclusions 

The conclusions for practical applications are quite 
straightforward: 

(1) In complex Langevin simulations one should use 
Ni = 0. If one wants to do a random walk simulation 
or an iterate of the Fokker-Planck operator, this is not 
possible, but one should make sure that Ni <C 1. 

(2) Always validate the simulations by comparison 
with other trusted computations for parameter values 
where other methods are available. 

(3) Wilson or Polyakov loops of higher charge will gen- 
erally have much larger fluctuations. In general they 
are not needed for physics applications, but they may 
be worth looking at because they can give information 
about the probability distribution. 
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Y 




(U) 




{V-) 


{U-") 


0.032 


CLE 


1 0.48331(047) 1 0.59265(058) 


0.12311(037) 1 0.19865(055) 


0.158 


CLE 


|0.48331(047) 1 0.59265(058) 


0.13211(037) |0.19865(055) 


0.474 


CLE 


|0.48331(047) 1 0.59265(058) 


0.13211(037) |0.19865(055) 


0.790 


CLE 


|0.48331(047) 1 0.59265(058) 


0.13211(037) |0.19865(055) 


1.581 


CLE 


1 0.48331(047) 1 0.59265(058) 


0.13211(037) 1 0.19865(055) 


CO 


CLE 


|0.48331(047) 1 0.59265(058) 


0.13211(037)|0.19865(055) 


exact 




1 0.48356 


0.59297 


0.13065 


0.19646 



TABLE I: Ni = 0.0 



Y 




(U) 


(IT-*) 


{V) 




0.032 


CLE 
FPE 


0.48344(046) 
0.48358 


0.59281(057) 
0.59297 


0.13154(036) 
0.13065 


0.19778(054) 
0.19646 


0.158 


CLE 
FPE 


0.48369(047) 
0.48361 


0.59313(057) 
0.59302 


0.13194(036) 
0.13065 


0.19835(055) 
0.19645 


0.474 


CLE 
FPE 


0.48370(047) 
0.48376 


0.59313(058) 
0.59307 


0.13188(036) 
0.13066 


0.19822(054) 
0.19639 


0.790 


CLE 
FPE 


0.48372(048) 
0.48405 


0.59319(059) 
0.59292 


0.13186(036) 
0.13076 


0.19812(055) 
0.19618 


1.581 


CLE 
FPE 


0.48382(047) 
0.48647 


0.59330(058) 
0.59033 


0.13183(037) 
0.13318 


0.19792(056) 
0.19449 


oo 


CLE 


0.48396(047) 


0.59345(058) 


0.13089(042) 


0.19685(056) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE II: iV> = 0.00001 



Y 




(U) 




{W) 




0.032 


CLE 
FPE 


0.48345(050) 
0.48366 


0.59280(062) 
0.59290 


0.13114(035) 
0.13066 


0.19714(052) 
0.19646 


0.158 


CLE 
FPE 


0.48385(054) 
0.48371 


0.59335(067) 
0.59313 


0.13153(039) 
0.13064 


0.19771(059) 
0.19644 


0.474 


CLE 
FPE 


0.48402(052) 
0.48399 


0.59357(064) 
0.59348 


0.13148(037) 
0.13059 


0.19762(056) 
0.19637 


0.790 


CLE 
FPE 


0.48425(053) 
0.48425 


0.59377(065) 
0.59379 


0.13146(038) 
0.13056 


0.19764(057) 
0.19627 


1.581 


CLE 
FPE 


0.48485(051) 
0.48476 


0.59454(063) 
0.59437 


0.13115(036) 
0.12363 


0.19725(054) 
0.19583 


00 


CLE 


0.48527(051) 


0.59504(063) 


0.12944(047) 


0.19387(079) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE III: iVj = 0.0001 



Y 




{U) 




{W) 




0.032 


CLE 
FPE 


0.48342(049) 
0.48373 


0.59293(060) 
0.59226 


0.13095(038) 
0.13074 


0.19692(057) 
0.19602 


0.158 


CLE 
FPE 


0.48410(051) 
0.48399 


0.59346(062) 
0.59344 


0.13068(037) 
0.13065 


0.19632(056) 
0.19647 


0.474 


CLE 
FPE 


0.48538(051) 
0.48488 


0.59497(064) 
0.59457 


0.13090(038) 
0.13050 


0.19702(058) 
0.19623 


0.790 


CLE 
FPE 


0.48619(051) 
0.48572 


0.59605(061) 
0.59559 


0.13066(041) 
0.13027 


0.19655(063) 
0.19589 


1.581 


CLE 
FPE 


0.48762(052) 
0.48729 


0.59761(062) 
0.59753 


0.12934(044) 
0.12934 


0.19503(064) 
0.19449 


00 


CLE 


0.48890(050) 


0.59942(063) 


0.12392(068) 


0.18602(086) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE IV: Ni = 0.001 



Y 




(U) 


([/"*) 


{W) 




0.032 


CLE 
FPE 


0.48047(049) 
0.48073 


0.58917(060) 
0.58856 


0.12935(035) 
0.12894 


0.19450(053) 
0.19328 


0.158 


CLE 
FPE 


0.48240(051) 
0.48313 


0.59150(063) 
0.59163 


0.13039(037) 
0.13055 


0.19613(056) 
0.19584 


0.474 


CLE 
FPE 


0.48747(053) 
0.48739 


0.59753(067) 
0.59760 


0.13080(041) 
0.13064 


0.19645(059) 
0.19646 


0.790 


CLE 
FPE 


0.49012(055) 
0.49020 


0.60078(069) 
0.60107 


0.12952(045) 
0.12984 


0.19437(069) 
0.19525 


1.581 


CLE 
FPE 


0.49532(052) 
0.49550 


0.60714(068) 
0.60758 


0.12647(054) 
0.12668 


0.19056(085) 
0.19051 


00 


CLE 


0.49995(049) 


0.61287(063) 


0.10733(202) 


0.16132(376) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE V: iVj = 0.01 



Y 




(U) 




{W) 




0.032 


CLE 
FPE 


0.46720(051) 
0.46763 


0.57290(063) 
0.57253 


0.12150(034) 
0.12118 


0.18269(051) 
0.18165 


0.158 


CLE 
FPE 


0.46985(054) 
0.47070 


0.57620(066) 
0.57629 


0.12335(037) 
0.12363 


0.18557(055) 
0.18533 


0.474 


CLE 
FPE 


0.48570(057) 
0.48608 


0.59533(071) 
0.59542 


0.13167(042) 
0.13184 


0.19766(065) 
0.19793 


0.790 


CLE 
FPE 


0.49664(055) 
0.49644 


0.60871(068) 
0.60858 


0.13166(046) 
0.13166 


0.19731(074) 
0.19799 


1.581 


CLE 
FPE 


0.51104(056) 
0.51055 


0.62580(080) 
0.62601 


0.12366(074) 
0.12363 


0.18370(114) 
0.18595 


00 


CLE 


0.52272(051) 


0.64013(069) 


0.07553(204) 


0.11359(334) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE VI: Ni = 0.05 



Y 




{U) 




{W) 




0.032 


CLE 
FPE 


0.45144(051) 
0.45203 


0.55358(062) 
0.55343 


0.11267(035) 
0.11236 


0.16941(053) 
0.16843 


0.158 


CLE 
FPE 


0.45456(055) 
0.45513 


0.55742(067) 
0.55722 


0.11453(038) 
0.11479 


0.17224(054) 
0.17207 


0.474 


CLE 
FPE 


0.47443(057) 
0.47501 


0.58151(071) 
0.58165 


0.12738(044) 
0.12759 


0.19117(065) 
0.19136 


0.790 


CLE 
FPE 


0.49560(053) 
0.49575 


0.60668(073) 
0.60740 


0.13448(052) 
0.13446 


0.20120(081) 
0.20198 


1.581 


CLE 
FPE 


0.52141(093) 
0.52157 


0.63883(093) 
0.63948 


0.12345(098) 
0.12426 


0.18569(144) 
0.18692 


00 


CLE 


0.54104(059) 


0.66230(082) 


0.05049(354) 


0.07723(522) 


exact 




0.48356 


0.59297 


0.13065 


0.19646 



TABLE VII: JV/ = 0.1 



